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ABSTRACT 

Measuring dark matter substructure within galaxy cluster haloes is a funda- 
mental probe of the ACDM model of structure formation. Gravitational lensing 
is a technique for measuring the total mass distribution which is independent of 
the nature of the gravitating matter, making it a vital tool for studying these 
dark-matter dominated objects. We present a new method for measuring weak 
gravitational lensing flexion fields, the gradients of the lensing shear field, to 
measure mass distributions on small angular scales. While previously published 
methods for measuring flexion focus on measuring derived properties of the lensed 
images, such as shapelet coefficients or surface brightness moments, our method 
instead fits a mass-sheet transformation invariant Analytic Image Model (AIM) to 
each galaxy image. This simple parametric model traces the distortion of lensed 
image isophotes and constrains the flexion fields. We test the AIM method us- 
ing simulated data images with realistic noise and a variety of unlensed image 
properties, and show that it successfully reproduces the input flexion fields. We 
also apply the AIM method for flexion measurement to Hubble Space Telescope 
observations of Abell 1689, and detect mass structure in the cluster using flexion 
measured with the AIM method. We also estimate the scatter in the measured 
flexion fields due to the unlensed shape of the background galaxies, and find 
values consistent with previous estimates. 

Subject headings: gravitational lensing, flexion, galaxy clusters, Abell 1689 
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Introduction 



Gravitational lensing is a well developed tool for measuring the mass distributions 
of galaxy clusters. Maps of the lensing convergence (the projected surface mass density 
between the observer and the source object scaled to the critical lensing density) in galaxy 
clusters, have been constructed using a variety of techniques for weak-field (or simply 



"weak") lensing observations (e.g., Kaiser et al. 1995 Hoekstra et al. 1998 Refregier & 



Bacon 2003 Kuijken 2006), multiple-image (or "strong") lensing observations (e.g., Kneib 



et al. 1993 Broadhurst et al.||2005 i Diego et al. 2005 Coe et al. 2010a), and combining both 



strong and weak lensing observations (e.g., Bradac et al. 2005). Joint weak+strong lensing 



approaches have the advantage of utilizing measurements of a broader range of lensing field 
strengths to produce a unified mass map. Lensing is a vital tool for measuring the total 
masses and mass distributions of galaxy clusters, as it directly probes the dark matter 
which dominates the mass budget of these large cosmological objects. 

While most work to date has focused on measuring image displacements and linear 
image distortions (strong lensing and weak lensing shear, respectively), gradients in 
the shear field which are significant on the scale of the lensed galaxy images produce 
second-order lensing distortions called "flexion." These lensing effects, which can be 
decomposed into a comatic distortion (1-flexion) and a trefoil distortion (3-flexion), are 
probes of small scale structures in the mass distribution. 

The two primary methods for measuring flexion in the current literature are shapelets 



(Goldberg & Bacon 2005 Massey et al. 2007; Leonard et al. 2007) and surface brightness 



moments/Higher Order Lensing Image Characteristics (HOLICs; Irwin & Shmakova 2006 



Okura et al. 2007 Schneider & Er 2008). The shapelet method is, in a certain sense, a 



refinement of moment methods, since shapelet coefficients are measured by calculating a 
specific combination of surface brightness moments (corresponding to Hermite polynomials) 
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using Gaussian weights. The inversion of shapelet coefficients (or simply "shapelets") 
to estimate the lensing fields utilizes a Taylor series expansion of the surface brightness 
profile to linearize the lensing transformation with respect to the shapelets basis and takes 
advantage of the well-described mathematics of the quantum simple harmonic oscillator. 
This linear approximation of the surface brightness profile is how the shapelets method 
fundamentally differs from moment methods, which use ratios of measured moments to 
approximate the flexion fields. In both approaches, assumptions on the unlensed surface 
brightness profile are made in order to measure the lensing fields. These assumptions are 
that odd-order shapelets/moments in observed images are sourced by lensing flexion alone, 
meaning that the intrinsic two-dimensional skewness of unlensed galaxies is assumed to 
be, on average, zero. The magnitude of the flexion scatter from intrinsic galaxy profile 
shapes is an important quantity to understand, as it sets the finite resolution of flexion 
mass measurements. This limit is analogous to how the distribution of intrinsic ellipticities 
limits the spatial resolution of weak lensing shear mass measurements. 

In each of the current flexion measurement methods, there are systematic issues 
that must be addressed. Shapelet-based lensing measurements have been shown to be 



unreliable and biased for large shears (Melchior et al. 2010). This is a significant concern 



for measuring flexion, since the shear can not be assumed to be small in the regime where 
flexion fields are measurable. Moment-based measurements of flexion rely on high-order 
surface brightness moments, making the strength of the lensing signal extracted very 
sensitive to the window functions used in the moment calculations: though the flexion 
symmetries under coordinate transformations are matched by third-order moments of the 
surface brightness, the normalization of the moments required to isolate the flexion signal 



includes both fourth- and sixth-order moments (Okura et al. 2008). Characterizing the 
noise properties of these moments and how the noise propagates into uncertainties in 
the measured flexion fields is a complex problem. Additionally, the implementations of 
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both these methods and their application to observational data have not addressed the 



well-established mass-sheet degeneracy ( Gorenstein et al. 1988 ; Saha 2000 ) in a consistent 



manner, as discussed by Schneider & Er (2008). These issues are not insurmountable, but 



further work needs to be done to refine flexion measurement, and alternative approaches for 
measuring flexion can provide a valuable insight to the strengths and weaknesses of each 
method. 

In this paper we introduce a new method for weak lensing flexion measurement. 
Instead of measuring derived quantities (such as weighted surface brightness moments, 
as in both the shapelets and HOLICs methods), we instead fit the lensed galaxy objects 
with a parameterized, Analytic Image Model (AIM) which is invariant to the mass-sheet 
degeneracy. By comparing model images to the data image in "pixel-space" and optimizing 
a figure of merit over a reasonable range of model parameter values, we constrain the flexion 
fields. This method has the advantage that surface brightness errors are well understood 
(and typically Gaussian), thus the optimization algorithm can provide reliable estimates 
of the errors on the best-fit parameters. Uncertainties in the flexion measured by shapelet 
and moment methods are quantified on average, rather than for each individual object, and 
primarily estimate the uncertainty in the mass reconstruction instead of the uncertainty 
in the measured flexion. Direct error estimates for each component of the 1-flexion and 
3-flexion values from each individual lensed object is a very desirable property, as it allows 
us to accurately weight the flexion measurements from each object in mass reconstructions. 

This paper is structured as follows: £j2] reviews the basic flexion formalism used 
throughout the paper. £j3] describes the the principle of the AIM method and the specific 
implementation used here. §|4] describes the procedure used to test the AIM method on 
simulated data images, validating the accuracy of the fitting procedure and the accuracy 
of the error estimates. §[5] describes the Hubble Space Telescope Advanced Camera for 
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Surveys Wide Field Camera (HST/ACS-WFC) observations of the galaxy cluster Abell 
1689, which is used as a proving ground for the AIM method in a real data application. We 
present detection of the substructure in A 1689 using only 1-flexion measured with the AIM 
method, and the structure agrees with previously measured lensing maps. In §[6] we discuss 
our estimate of the scatter in flexion measurements due to the intrinsic (unlensed) shape of 
galaxies. 



2. Flexion Formalism 

The distortion of images by gravitational lensing is described by a coordinate 
transformation between positions 9 in the image plane and positions (5 in the unlensed 
source plane 

p = e-vm, (i) 

with 

4(0) = ^ J d 2 0' K(0')\n\0-0'\. (2) 

We use complex notation here, where coordinates and lensing fields are denoted by complex 
quantities x = X\ + ix%, and the derivative operator, V = d\ + 182-, is also complex, k and ijj 
are real-valued, but all other lensing quantities have non-trivial imaginary components. 



2.1. Second-Order Local Lensing 

In the case of weak lensing, the transformation in Equation [T] is expanded in a Taylor 
series valid in the neighborhood of the position of the image (#0 in the image plane, Po in 
the source plane) in question. To include the flexion fields, this expansion must extend to 



-7- 



quadratic order, and yields 

0o + & = 9 - Wk + 

(l-/c)0- 7 0*- (3) 

-TV - -Fee* - ^g(9*) 2 . 

The first-order lensing fields (convergence and shear) are k = and 7 = \V 2 ip; 

the-second order fields (flexion) are J 7 = Vk = V*7 and Q = V7. Throughout this paper 
we will refer to J 7 as 1-flexion and Q as 3-flexion, refering to the spin symmetry of the fields 
with respect to coordinate rotations. Elsewhere in literature, these lensing fields are called 
"first flexion" and "second flexion" , respectively. We prefer the spin-based notation because 
it indicates a physical property of the flexion fields, rather than an arbitrary ordering. 
All derivatives of the lensing potential are evaluated at 9q. This notation follows that of 



Schneider & Er (2008), which simplifies the tensor notation of Goldberg & Bacon (2005). 



When considering a single lensed image, the constant terms in Equation [3] are 
degenerate and unmeasureable, and can be neglected in favor of small deviations 9 about 
the observed image center 9 . This gives a second-order local lensing equation 

=(1 - K )e - 7 0*- 

111 W 

-t*9 2 - -Tee* - -Q{9*) 2 . 

This local lensing equation is sufficient for producing "arced" images from intrinsically 
circular or elliptical ones. It is valid in the regime where the dimensionless products of the 
image size and the flexion fields, ai\J : \ and are small with respect to unity. The 

effect of flexion becomes significant when ai\J- r \ and/or aj\Q\ are small but non-negligible 
corrections to the linearized lensing equation. If the lensing flexion is significant, the weak 
lensing approximations, k <C 1, I7I <C 1, and (3(9) is linear, are typically no longer valid. 
This has significant implications for characterizing the lensing transformation in regards to 
the mass-sheet degeneracy. 
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2.2. Mass- Sheet Invariance 



The mass-sheet degeneracy is present in all local lensing measurements (Gorenstein 



et al. 


1988 


Saha 


2000) 



of lensed images at a single redshift can only be constrained to a family of solutions 
{k', 7', J 7 ', Q'}, where k' = Xk + (1 — A), 7' = A7, T' = AJ 7 , and Q' = \Q are equally 
valid solutions for any non-trivial value of A. This degeneracy corresponds physically to an 
unconstrained source plane scale. Multiply-imaged (strongly lensed) sources observed at 
multiple redshifts are required to constrain the source plane for those objects, fixing the 



constant terms in Equation 3 and breaking the mass-sheet degeneracy (Falco et al. 1985 



Schneider & Seitz 1995). 



There are combinations of the lensing fields which are invariant to mass-sheet 
transformations. These are the reduced shear g, reduced 1-flexion and reduced 3-flexion 
\l/3. By referring to a rescaled, fiducial source plane where /3' = 0/(1 — k), rather than the 
true source plane, the lensing transformation becomes 



= 9-g6* - -*^ 2 - ^#100* - -# 3 {tf :i ) 



*\2 



(5) 



where g — 7/(1 — k), ^1 = — k), and ^3 = Qji\ — k). This notation is similar to 



that of Schneider & Er (2008), though we do not absorb a factor of 1/4 into the definition 



of ^1 and \& 3 to better match our parameter values to the flexion estimators of Er et al 



(2010). This formulation quantifies the lensing transformation in terms of three complex 
variables rather than three complex and one real, casting the lensing transformation in 
terms of mass-sheet transformation invariant fields only. Each of the three reduced lensing 
fields is specifically defined and measurable for single images, though determining k from 
these reduced lensing fields requires a constraint of the mass-sheet degeneracy. From this 
point on we will use the transformation in Equation [5] and drop the prime notation from f3' 
simplicity. 
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3. Measuring Flexion with AIM 

In this section we describe our alternative to shapelets and moments for measuring weak 
lensing flexion. Rather than measuring either set of derived quantities from the observed 
images and making assumptions on their properties in the (unobservable) unlensed image, 
we instead forward-model the observed image analytically. By assuming a parameterized 
ansatz for the unlensed image and applying the quadratic lensing equation in Equation [5j 
a model image is generated and compared to the observed, lensed image. As we show, the 
ansatz need not be exactly matched to the true galaxy shape: an elliptical Gaussian is 
sufficient for most objects. This forward-modeling approach is what sets the AIM method 
apart from shaplet and moment based methods, which instead attempt an inversion of 
image characteristics to measure 1- and 3-flexion. 



3.1. The Analytic Image Model Method 

As a consequence of Liouville's theorem, the surface brightness of an astrophysical 
source is conserved by gravitational lensing. This means that the observed surface brightness 
Jobs at an image-plane position 9 in terms of the source plane surface brightness J src is 

/obs [9] = Jsrc W)\ , (6) 

where (3(9) is the lensing coordinate transformation. If we assume that the intrinsic, 
unlensed surface brightness profile can be well-described by a set of model parameters 
{Pint}, then the lensed model image is defined to be 

^AIM [9] Pint, Plcns] = -^src [/3(6, Plena), Pint] , (7) 

where pi ens is a set of parameters characterizing the lensing transformation. For 
measuring flexion with such an analytic image model, the parameter set pi ens = 
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{gi, g 2 , ^11, ^12, ^31, ^32} defines this transformation. Here gi, ^n, and ^31 are the real 
parts, and g 2 , ^12, and \& 32 are the imaginary parts of g, ^1, and ^3, respectively. 

The AIM is optimized by minimizing the figure of merit 

2( , V- (/ob S (^ (n) ) - /AIM(^ (rl) ;Pin t ,Plens)) 2 , Q , 

X faint, Plena) = 2^ 2 ' ( 8 ) 

n ° n 

where 9^ is the image-plane position of the n th pixel and a\ is an estimate of the variance 
in that pixel's value. The parameter set which minimizes this figure of merit is our estimate 
of the true set of intrinsic and lensing parameters. 



3.2. Implementation 

We implement the AIM method using a Levenberg-Marquardt minimization algorithm 
called MPFIT ( |Markward"t1|2009[ ). MPFIT is written in the Interactive Data Language^ (IDL). 



Each model parameter is restricted to a fixed range of values encompassing the plausible 
span for that parameter. The range for each parameter is described below and listed in 
Table ffl 

The Levenberg-Marquardt minimization algorithm yields not only best-fit parameters, 
but also a full covariance matrix (Cjk) for all parameter pairs. Both the variance estimate 
for each parameter (crj = Cjj) and the correlation matrix (pjk = Cjk/ \jCjjCkk) provide 
additional indicators for the quality of the fit. The variance estimates directly measure 
the constraint that the AIM minimization places on each model parameter for each object 
analyzed, and the correlation matrices can be used to evaluate and improve the model 
parametrization we use. 

Models with many parameters (this implementation has twelve) can often have 



1 IDL is produced by ITT Visual Information Solutions: 



http:/ /www. ittvis.com 
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significant parameter degeneracies, making the correlation matrix obtained from the 
fitting a valuable tool for evaluating the specific model parametrization. Large correlation 
coefficients (\pjk\ ~ 1) indicate significant degeneracies which can prevent the minimization 
algorithm from efficiently and accurately converging. 

There is a lensing-parameter /shape-parameter degeneracy which is well known to 
the lensing community: the shear/ellipticity degeneracy utilized in standard weak lensing 
studies to measure shear. Our approach to mitigating this degeneracy, fixing the shear 



parameters during fitting, will be discussed in £3.3, and the simulations described in §[4] 
confirm that there are no other strong degeneracies. Modest correlations \pjk\ ~ 0.7 exist for 
some lensed images and parameter combinations, but even in these instances convergence 
is robust and accurate. 

The lensing transformation is characterized by six variables: the three complex, 
reduced lensing fields g, and ^3 from Equation [5} The range allowed for each these 
parameters is listed in Table [TJ We fix the shear during fitting to address the degeneracy 
between the shear and the intrinsic ellipticity, and so do not set explicit limits on the range 
of values for the two shear parameters. §3.3| discusses this point in more detail. 



Though it is well-known that intrinsic galaxy profiles are not Gaussian (see, e.g., 
Graham Driver||2005 and references therein, for a review of non-Gaussian galaxy shapes), 



we assume an elliptical Gaussian model for the unlensed images. This is jusitified because 
the goal of the AIM method is not to match the profile of the galaxy image exactly. Rather, 
the goal is to model the distortion by lensing of the galaxy isophotes. A Gaussian profile 
effectively acts as a simple window function that identifies the isophotes with a minimal 
set of parameters. The assumption that the unlensed isophotes are elliptical is a standard 
lensing assumption. For an elliptical unlensed image, the isophotal ellipses will be similarly 
distorted for a variety of intrinsic surface brightness profile slopes. Simulations which 
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assumed Sersic profiles for the unlensed image model failed to converge at a much higher 
rate ~ 25% and yielded flexion estimates with errors nearly an order of magnitude greater. 
This is due to degeneracies between the normalization, size and Sersic index parameters. 

With an elliptical Gaussian ansatz for the unlensed image, there are twelve model 
parameters: six for the unlensed profile and six for the lensing transformation. The 
Gaussian parameters are combined into two real-valued variables, log So and a; and two 
complex variables, 6 C and e. In this parametrization, the unlensed surface brightness profile 
is defined by 

with 

r(/3) 2 = (1 + ee*)/3/3* - /3 2 e* - (/3*) 2 e. (10) 

The image scale and complex ellipticity e are defined in terms of the semi-axes and position 
angle (A, B, and 0) as 

a — VAB; e=^|e 2 ^ (A > B). (11) 
The center position of the image is included by defining 

/AIM [0] = /src W ~ e c )) . (12) 

We define the center position in the image plane for a better a priori limit on the allowed 
parameter range for the center position, and the origin of the source plane coordinates is 
defined to be at the center of the unlensed image. Note that this center position does not 
refer to the image centroid, but rather the location in the image plane where the peak 
surface brightness of the source plane is mapped. Table [T] lists the parameter ranges allowed 
during minimization, which are chosen to be very broad and thus prevent parameters 
from reaching these limits. However, if a parameter does happen to reach the limits while 
optimizing, that fit is flagged as a non-convergence. 
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In general, there will be a non-negligible background surface brightness and a finite 
point-spread function (PSF) size. Since these observational effects are not well constrained 
by the individual lensed galaxy images, we take them as given for the AIM optimization. 
In practice, they will be determined externally to the AIM fitting; see ^5] for our approach. 



3.3. Shear/Ellipticity Degeneracy 

As is well known from weak lensing studies, there is a complete degeneracy between 
shear and ellipticity. Assuming only a linear lensing equation (neglecting flexion), the 



observed ellipticity is given by (Bartelmann & Schneider 2001): 



e °bs = ^— " for \g\ < 1, 

1 + eg* 

(13) 

= for \g\ > 1. 

e* + g* 

Except in the limit of extreme image distortions, this degeneracy prevents the shear from 
being well-constrained by a single lensed image. This is why weak lensing shear studies 
must average the ellipticities of nearby galaxies in order to infer the lensing shear. For 
the AIM measurements of flexion presented here, the presence of this degeneracy means 
that the four-dimensional parameter volume spanned by e and g only has two constraints 
(e Q bs)- In order to accurately constrain flexion using the y 2 minimization, we must reduce 
this part of the AIM parameter space by two dimensions to prevent the minimization 
algorithm from converging to degenerate local minima rather than finding the true global 
minimum. We do not expect that the specific constraint on e and g that we choose will 
affect the flexion estimation, since the flexion fields distort images with different symmetries 
than shear and ellipticity (spin-1 and spin-3 flexion, versus spin-2 shear and ellipticity); 
and indeed, the moderate correlation coefficients between shear /ellipticity parameters and 
flexion parameters in our simulations (described in Q bear out this expectation. 
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To remove this degeneracy we choose to fix the two shear parameters to an assumed 
input model. This eliminates two dimensions from the parameter space, drastically reduces 
the convergence time of the fitting algorithm (by a factor of ~ 10), and allows the flexion to 
be accurately and precisely constrained. In simulations where both the ellipticity and the 
shear are allowed to vary freely, the flexion parameter values returned by the minimization 
place no significant constraint on the lensing flexion, and nearly all the parameter values 
are consistent with zero flexion within the error estimates derived from the parameter 
covariance matrix. 

We choose to fix the shear, rather than fixing the ellipticity, because it is possible 
to select a specific a priori shear model which is observationally motivated by additional 
data on the galaxy cluster in question. The shear can be estimated from mass proxies 
(e.g. standard weak lensing shear analysis, X-ray temperature or luminosity, or optical 
richness) and the assumption of a standard mass profile (e.g., a non-singular isothermal 
sphere or NFW profile). Significant departures from the input shear model would be 
measurable as departures from the expected distribution of ellipticity magnitudes and 
orientations in the lensed galaxies. 

Flexion measurements are a natural addition to joint weak+strong lensing mass 
mapping methods, as flexion directly constrains the mass gradients. Mass reconstruction 



formalisms, such as that of Bradac et al. (2005) or the non-grid-based method of Deb et al. 



(2008), can be easily modified to utilize data which further constrains the lensing potential. 
As a less-sophisticated method which does not use direct shear estimates, an iterative 
process of integrating flexion measurements could be used to correct the assumed input 
shear model (and thus the mass model) for the effects of lens substructure. The resulting 
updated shear model could then be used as an input to the next iteration of flexion fitting, 
and the process repeated until the shear model converges. This type of procedure would be 
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a test of the self-consistency of the measured shear and flexion fields; however it is beyond 
the scope of this work. The simulations we present in §|4] confirm that our basic approach to 
the shear /ellipticity degeneracy (fixing the shear and fitting for the ellipticity) is adequate 
for successfully measuring the lensing flexion, even when significant errors are present in 
the assumed shear field. 

4. Testing the AIM Method 

We test the efficacy of the AIM method for measuring flexion using simulated data 
images. Image noise and pixel scale are chosen to match the Abell 1689 data described 
in £|5j Data images are generated by randomly selecting intrinsic shape and lensing 
parameter values, assuming a specific intrinsic galaxy profile and applying the quadratic 
lensing transformation. These images are circularly windowed such that the observed 
image centroid is situated in the center of the simulated data image, to emulate images 
that would be extracted from real data. The AIM method is used to determine best-fit 
parameter values. The lensing parameters are either selected uniformly from fixed ranges 
or from a Singular Isothermal Sphere (SIS) lens model constrained to the same volume of 
lensing parameter space. A set of 1000 simulated images are generated to explore the AIM 
parameter space for each lensing parameter selection method. 

4.1. Intrinsic Shape Parameter Selection 

Each intrinsic image follows a Sersic profile, 

I oc exp(-0.5(r/a) 1/n ), (14) 

with the index n uniformly selected from the range 0.2 to 5. For reference, n = 0.5 is a 
Gaussian, n — 1 is exponential, and n = 4 is a de Vaucouleurs' profile. The intrinsic size 
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parameter a is drawn from a uniform range between (X'25 and 1'.'25. The normalization 
log So (in arbitrary units) is selected uniformly from a range such that the object would be 
a 2-<T detection or better. 

The ellipticity magnitude is drawn from a Gaussian distribution of the form 



I f-\e 

72 



P(e) = — 27i / ; 7 2V\ exp — 2~ ' ( 15 ) 

7ra z (1 — exp(— 1/er J) \ 



with a = 0.2, but truncating the distribution at |e| = 0.9, as in Schneider (1996). The phase 
of the ellipticity is uniformly selected to be from to 2tt, randomly orienting the unlensed 
ellipse. The center position 9 C is selected such that the image centroid is located at the 
center of the data image. 

The AIM method assumes that the intrinsic profile is Gaussian, while the simulated 
data images explicitly include non-Gaussian Sersic profiles. As a result, the best-fit values 
for log Sq and a will not match well to the input values if the Sersic index differs significantly 
from 0.5. However, the shape of the isophotes is still constrained and the other intrinsic 
shape parameters (ellipticity and center position) can be well fit. The exact true values of 
log So and a for the unlensed image need not be determined for an accurate lensing analysis. 



4.2. Lensing Parameter Selection 

Lensing parameters are randomly selected using one of two methods. The first method 
is to simply draw each parameter from a uniform distribution in a volume of parameter 
space. For this uniform selection, we restrict the parameters to \g n \ < 0.5, and \^f mn \ < 0.1 
arcsec -1 . Images are generated using the quadratic lens equation. 

The second selection method involves using the uniformly selected parameters to 
choose a specific SIS lens model. We do this to test the efficacy of the AIM method on 
images which have been lensed analytically from a well-defined lensing potential, rather 
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than restricting the lensing to be strictly quadratic. For an image located at a radius r and 
position angle (f) from the center of a SIS, the lensing fields are 

^ = -^Vexp(2z0), (16) 
zr — ue 

*i = - J E ft ^ exp(^), (17) 
r(2r — #e) 

and 

* 3 = K^) exp(3 ^ (18) 

where is the Einstein radius of the SIS. From the uniformly selected lensing parameters, 
we calculate t/6e from \g\\ <fr from \l/ 12 /\l/ 11 ; and r and 9e are separately determined 
by combining with \g\. This allows us to select a uniform distribution of lensing 
parameters while requiring that the lensing field combinations be those of a SIS. 



4.3. Image Generation 

Simulated images are created and windowed to include only those pixels at a distance 
from the observed image centroid of a factor of 1.5 times the observed semi-major axis 
length A. This window size is empirically determined to be sufficiently large to measure 
flexion accurately while being independent of ellipticity and not being so large as to include 
excessive sky noise. This choice of 1.5 A is motivated by similar concerns presented by 



Goldberg & Leonard (2007). For uniformly selected lens parameter images the simulated 



data images are generated using the AIM described in §3.1| For SIS lens parameter images, 
the lensing transformation is known exactly, and so the images are analytically lensed 
without approximating the transformation to second order. 

Each data image is convolved with a circular, Gaussian PSF with a size similar to 
that of the HST/ACS-WFC. Gaussian-distributed pixel noise is included with a standard 
deviation chosen to match the dataset described in $5l We do not use a more accurate HST 
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PSF model (e.g. TinyTim; Krist 1993) because the effects of the PSF on the measured 
flexion are expected to be negligible due of the small PSF size (0'.'08 FWHM) relative to the 
size of galaxies which we select for analysis (> (X'24 FWHM). Assuming that the PSF has a 
shape corresponding to a 1-flexion J-psF, the effect on the measured 1-flexion scales as 

•^induced ~ ^~PSF^ ~~4 (19) 

a PSF ~i~ a image 

sharply damping out PSF-induced flexion for galaxies only modestly larger than the PSF 



(Leonard et al. 2007). PSF ellipticity couples to the ellipticity and shear parameters, not to 
the flexion parameters. The approximation of a Gaussian PSF reduces the computational 
intensity of the AIM method and is a valid simplification, though a more complex PSF 
model would need to be used for images where the PSF size were larger, such as for 
ground-based imaging. 



4.4. Initial Parameter Selection 

Initial parameter values for intrinsic shape parameters are determined from surface 
brightness moments of the data images. The center position is initially set at the observed 
centroid of the simulated data image. The flexion parameters are all initialized at zero, 
meaning that the "first guess" model image is always an elliptical image, and non-zero 
flexion measurements will be positive detections by the fitting algorithm. 



As noted in £3.3, there is a degeneracy between the ellipticity and the shear, requiring 
that the shear parameters be fixed to some set value for the minimization to converge. 
For these simulations, each shear parameter is fixed to a value deviated about the true 
value by a Gaussian distribution with a standard deviation of 0.2 (no units). Accordingly, 
the initial ellipticity parameter values are adjusted using Equation [13] so that e Q b s matches 



the ellipticity determined from the simulated data image moments. This ensures that the 
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observed ellipticity of the initial model image matches that of the data image. 

4.5. Simulation Results 

The AIM method converges to a minimum rapidly and consistently (> 99% 
convergence, typically in less than 100 iterations). The final figure of merit \ 2 P er Degrees 
of Freedom (D.o.F.) for converged fits have mean values of ~ 1.0-1.2. A typical data image 
has 500-1000 D.o.F. Successful, converged fits are identified from the larger ensemble as 
those having x 2 /D.o.F. < 1.5 and a(^ mn ) > 0.001 arcsec -1 . The latter condition is slightly 
couterintuitive, but it is useful for excluding fits which tightly constrain one or more flexion 
parameters to erroneous values in local minima of the x 2 surface. We are able to identify 
this threshold because we know the true input flexion values. Fits with one or more flexion 
parameter constrained this tightly by the minimization do not return accurate flexion 
estimates. 

The error between the true and fit normalization and size parameters (log So and a) 
are small, with RMS values of 0.1 for logS an d / .'04 for a. These parameters are positively 
correlated, with typical correlation coefficients of p\ og s /a ~ 0.7. For low surface brightness 
Sersic input profiles with indices significantly different from the Gaussian value of 0.5, 
the errors in logS an d a become exaggerated. This scatter from model mismatching 
is expected, and residual images show that the central region of the image is not well 
estimated. The flexion can still be accurately determined, as flexion primarily distorts 
the isophotes away from the image center. This property is indicated by the absence of 
any significant correlations between the flexion parameters and either the normalization or 
size parameters. The RMS center position error is ~ 0'.'05 (approximately one pixel). The 
ellipticity is measured with an RMS error of ~ 0.2, matching the distribution of the input 
shear values. Similar simulations fixing the shear parameters to their true values yield an 
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RMS error on the ellipticity of ~ 0.01. Figure [I] shows a comparison of fit and input 1- and 



3-flexion values from simulations. As has been observed in earlier work (e.g. Goldberg & 



Leonard 2007), there is more noise and scatter in measurements of 3-flexion. 



Flexion parameters are also well-recovered by the AIM method. Though there is some 
correlation between the flexion and ellipticity parameters, it is modest {\pjk\ ^ 0.7) and 
does not preclude accurate recovery of the input flexion fields. Typical error estimates from 
the flexion variances yield a(^ mn ) ~ 0.01 arcsec -1 , and the RMS deviation from the input 
flexion values is 0.007 arcsec -1 . This indicates that for those fits selected from within the 
acceptable ranges of x 2 /D.o.Y. and cr(\l/ mn ), the parameter variances provide an accurate 
estimate of the deviations of the measured flexion values from the true lensing fields. 



5. Application to A1689 



Flexion has been used to measure mass structure in two contexts to date: galaxy-galaxy 



flexion in field galaxies (Goldberg & Bacon 2005 Velander et al. 2011) and in the galaxy 



cluster Abell 1689. This massive cluster is one of the best-studied gravitational lenses, 
making it an ideal testbed for a new lensing analysis technique. Additionally, it is the 
only galaxy cluster to date where flexion has been used to measure the mass distribution. 



Leonard et al. (2007) used shapelets to measure flexion in A1689 and constrained mass 



structures with a parametric model. Okura et al. (2008) used HOLICs to measure flexion 



and found similar structure using weak lensing shear and flexion measurements and a Fourier 
mass reconstruction. In this section we apply the AIM method of flexion measurement to 
substructure detection in A 1689 as a real- world validation of our method. 

It is important to note that in both previous flexion analyses of A1689, the authors 



mistreat the mass-sheet degeneracy in some way. As shown by Schneider & Er (2008), this 
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is not a valid approximation for the interior region of the cluster where k ^ 1. This means 
that direct, quantitative comparisons between this paper and previous work are not feasible. 
Instead, we look for qualitative morphological similarities between the structures detected, 
with the intent of performing a direct comparison of flexion measurement methods in future 
work. 



5.1. The Data 

The data were obtained from the Hubble Legacy Archive^] (HLA) and reduced 
using the standard HLA pipeline. These data were originally obtained for proposal 
HST-GTO / ACS9289 (PI Ford). The data include deep observations in F475W (9.6 ks), 
F625W (9.6 ks), F775W (11.8 ks), and F850LP (16.6 ks). The observed field is 3^4 x 3(4, 
and 1" on the sky corresponds to 3.1 kpc at the redshift z = 0.187 of the cluster, assuming 
a standard cosmology. 



5.2. Object Selection 



Candidate objects for flexion analysis are selected using SExtractor (Bertin & Arnouts 



1996) in a two-pass strategy similar to those of Rix et al. (2004) and Leonard et al. (2007), 



with some modifications. The first pass is intended to select the large, bright, cluster 
member galaxies and foreground objects. These objects are then removed from the image in 
a process we describe below. The resulting "cleaned" image is then used in the second pass 
to identify the smaller, fainter, lensed background galaxy images. In all cases, a co-added, 
four-filter image is used for object detection to increase the signal-to-noise ratio of faint 



2 The Hubble Legacy Archive is located at: 



http: / /hla.stsci.edu/ 
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object detections, and the individual filter images are used for photometry. SExtractor 
uses the input data images and pixel weights to produce a background image, noise image, 
and object image for each of the four filters. The background and noise images contain 
position-variable estimates of the background surface brightness and the total surface 
brightness variance across the field. The object image is the same as the input image except 
that any pixels not determined to be associated with a detected object are set to zero. 

In the first pass, we select objects of large area (area > 50 pixel 2 ) to a low threshold 
(SNR > 1). This identifies object-associated pixels down to the background level, including 
the faint Intra-Cluster Light (ICL) associated with the central cluster galaxies. The detected 
objects are matched with known foreground and cluster member objects from literature 
( |Duc et al.|[2002[ |Coe et al.||2010b[ and the NASA/IPAC Extragalactic Database Q. All 



known objects at a spectroscopic redshift z < 0.2 are selected to be cleaned from the image. 
All pixels associated with any known object location, as determined in the SExtractor 
object image using a friends-of-friends algorithm, are replaced with Gaussian-distributed 
noise. The object pixels are cleaned in the four-filter detection image and in each individual 
filter image as well. The standard deviation of the noise, as determined by SExtractor and 
output to the noise image, is enhanced in the object pixels by a factor of 2. This is to 
account for statistical error in the subtraction of the image pixels as well as any systematic 
error in the subtraction process. 

The cleaned images are then re-run through SExtractor, though now allowing for 
slightly smaller objects (area > 25 pixel 2 ) with a more exclusive noise threshold (SNR > 3). 
Note that a new noise image for use in the flexion analysis is created in the second pass. 
This returns a catalog of objects for flexion analysis, along with initial parameter values for 
log So, a, and e. Those objects with a < 2 pixels are removed from the catalog as they are 



3 NASA/IPAC Extragalactic Database (NED): http://nedwww.ipac.caltech.edu/ 
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either spurious detections of noise peaks, stars, or are PSF-dominated galaxies which will 
have unreliable flexion estimates. Our flexion analysis catalog includes 764 objects in the 
field (a density of ~ 66 arcmin -2 ). 

There is an additional surface brightness cut which is applied after the fitting, removing 
objects with a peak surface brightness of S /(27ra 2 ) < 5 e~/s /square-pixel. This removes 
any remaining spurious noise detections and objects with low signal more effectively than 
using the SExtractor determined initial values. Application of this cut before the flexion 
analysis step removes ~ 30% of the usable, flexion-yielding objects in addition to the 
unsuitable objects. This is because for detections which are noise peaks or excessively 
faint, the AIM fitting algorithm naturally reduces log So and increases a so that the surface 
brightness falls below this threshold. For usable objects, there is a floor for log So and a 
ceiling for a beyond which the x 2 values become unfavorable, and so these faint objects 
remain in our sample and can be used to constrain the flexion fields. 



5.3. Flexion Fitting 



"Postage-stamp" images of each of the objects selected for flexion analysis are excised 
from the F775W data image and windowed to within a circular radius of 1.5 times the 
SExtractor-determined semi-major axis size about the observed centroid position. As 



noted in £3.3, the input shear parameters must be fixed for the fit to converge accurately. 



We choose a non-singular isothermal sphere (NIS) as a model for the cluster shear. The 
(non-reduced) shear and convergence for a NIS are 

E 



k{9) 



+ 9 2 



c 



and 



7(0) = k(9)9 2 



29r 



91 - 99* - 291 



*\2 



(20) 



(21) 
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We set 6e =49 / .'5 and 9c =17", as measured by Broadhurst et al. (2005). As noted by 
Broadhurst et al., it is unclear whether a NIS model, a NFW model, or some other mass 
model best describes the large scale cluster potential, particularly given the amount of 
known substructure in the cluster, but the NIS model is consistent with previous data and 
is a simple model to use for the shear. Since our approach focuses on measuring flexion 
rather than accurately modeling the shear, even with errors the assumed shear model is 
an adequate approximation, though certainly not a perfect representation of the true shear 



field. As noted in £3.3, the ellipticity will compensate for errors in the shear model. An 
aphysical model with the shear assumed to be uniformly zero across the entire field was 
also tested, with very few significantly different flexion measurements, though the median 
X 2 /D.o.F. was slightly elevated. This supports our assertion that the flexion parameters are 
minimally correlated with either the shear or ellipticity parameters. 



As in £j4| we remove all fits with x 2 /D.o.F.> 1.5, and those with (r($ mn ) < 0.001 
arcsec -1 . Figure [2] shows three example fits: one well below the x 2 /D.o.F. threshold, one 
just accepted, and one which was rejected. The typical rejected objects are either faint 
objects located near a significantly brighter one, or extremely low surface brightness objects. 
A final flexion catalog of 301 objects remained after these quality cuts were imposed. We 
use these 301 objects to detect the mass structure in A1689. For reference, Table [3] lists the 
position and fit parameters for the 50 objects with the largest 1-flexion signal-to-noise ratio. 



5.4. Mass Structure Reconstruction 



We reconstruct the mass structure from the measured 1-flexion field, using a modified 



version of the mass reconstruction technique introduced by Leonard et al. (2010). Their 



method, which is similar to the aperture mass measurement techniques used for shear (e.g. 



Schneider 1996), relates the weighted integral of the convergence on a circular aperture of 
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radius R about a center position 9 t 



o 



M ap (6 ) = / d 2 k(0 + o M|0|), (22) 

./[0|<fl 

to an integral of the 1-flexion: 

M ap (e )= [ d 2 e F E {9-e Q )Q T {\e\). (23) 

J\e\<R 



The E-mode 1-flexion, 



is the component of the 1-flexion which is radially oriented from 9 towards 9q. The 1-flexion 
kernel Qj is related to the aperture weight function w by the relations 

i r 

Q F ( X ) = — / w(x') x'dx', (25) 
x 

and 



o 



w(x) = - l -Q T {x) - (26) 

The kernel Qj has units of angle, T has units of inverse angle, while M ap and k are unitless. 
The weight function is constrained to go to zero smoothly at the aperture boundary and 
also to have the property that 

w(x') x'dx' = 0. (27) 







This latter restriction fixes the constant of integration to zero, eliminating the degeneracy 
k — y k -\- k,q from the mass reconstructions, where kq is a constant offset. Note that this 



is not the mass-sheet degeneracy. Leonard et al. (2010) assume that the flexion fields 



measured are the non- reduced flexion fields T and Q. 

However, as noted above, we can not measure J 7 , and instead measure The 
relations between M ap , k, and J 7 , as well as those between w and Qjr derive from the 
gradient relation between k and T = Vk. has its own gradient relation: = VK, 
where 

K = — In |1 - k\. (28) 



-26- 



We define an aperture statistic similar to M ap using the E-mode reduced 1-flexion, 
^ie, denned analogously to Te- This statistic is 



K, 



\<R 



\6\<R 



d 2 e K(9 + e )w(\e\) 

d 2 9 ^ 1E (9;9 )Q^ 1 (\9\), 



(29) 



where is related to w in the same way that Qjr is. Neglecting the difference between 
J 7 and (and thus the difference between M ap and K ap ) is a likely contributor to the 



anomalous B-mode convergence measured by Leonard et al. (2010). We use K ap , calculated 



from our measured mass-sheet invariant 1-flexion, to detect mass structure without making 
any assumptions about the redshift distribution of the source objects. 

We reconstruct K ap (9) from our flexion measurements using the following method. The 



integral in Equation 29 can be converted to a discrete sum over N flexion measurements 
simply: 



K ap (9 ) = J2^^(9^;9 )Q^(\9^\)q n . 



(30) 



Here q n is a normalized weight factor given by the inverse square of the 1-flexion errors 
returned by the AIM fitting. 



For the aperture weight function, we use the polynomial weights described in Schneider 



et al. (1998) and Leonard et al. (2010). These form a non-optimized, general family of 



polynomial functions characterized by a polynomial index I and a maximum aperture radius 
R. They are defined by the 1-flexion kernel 

2 + 1 



Q(r) = -At 



2tt 



r 1 



r 2 \ 1+/ 



and the aperture weight function 



w(r) = Ai 



(2 + If 



7T 



R 2 ) 



m) \2--i 



(31) 



r 

R 2 



(32) 
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with the normalization factor 



At 



4 T(7/2 + /) 

A r(3 + /) 



(33) 



chosen so that 




Q(r)2irrdr = 1 



(34) 



Because the aperture radius R can be freely selected, and the width of the kernel varied 
by changing I, this family of functions is sensitive to a wide range of mass structure scales. 
Features which can be identified in several mass reconstructions with different aperture 
radii and different slopes can be associated with real structure and not noise peaks. 

We use the apertures R =45", 60", 75", and 90", and polynomial indices I =3, 5, 
and 7. Figure [3] shows plots of w(r) and Qq ll (r)/R versus r/R (so scaled to make them 
unitless quantities). Higher polynomial indices make the kernel sensitive to a smaller range 
of radii within the aperture. This implies that high-index kernels will detect smaller mass 
structures, though will also be more susceptible to noise fluctuations due to the finite 
number of flexion samplings within the aperture. Lower indices reduce the amount of 
noise from discrete sampling, but also have a coarser resolution for mass structure and will 
therefore smooth out smaller structures. As the aperture radius increases, the flexion signal 
correlated with the center of the aperture falls off quickly, again smoothing over small 
structures while large structures continue to be detected. 

Because the systematics of this aperture statistic need to be further investigated, we 
evaluate the significance of the structure detected by constructing signal-to-noise ratio 
(SNR) maps from our flexion measurements instead of simply examining the K ap maps. 
This is analogous to the use of "B-mode" convergence maps in other weak-lensing studies 
to evaluate the significance of mass structure detections, though it is more directly related 
to specific parameter uncertainty estimates. 

The procedure for producing the SNR maps is as follows: For each combination 
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of R and I, we first construct a map of K ap on a 500x500 grid of points 6*0. We then 
construct 1000 "deviate" maps. For each deviate map, both components of the 1-flexion 
measurements for each object are shifted by a randomly-selected, normally-distributed 
value with standard deviation equal to the AIM fitting error for that flexion parameter. For 
each grid point 9q, the standard deviation of the deviate K ap values at that grid point is 
used to calculate the SNR map value 

SNR(^o) = ^4^r. (35) 

This yields maps of mass structure for each combination of R and I. Figure [4] shows 
SNR contours overlaid on the F425W/F625W/F775W color image of A1689 for R = 60" 
and I =3, 5, and 7. The major structures observable here are also detected in the other 
aperture sizes, though the larger apertures smooth over the smaller structures and reduce 
the magnitude of the SNR peaks. 

The SNR maps do not represent quantitative statistical significance values, as we expect 
that there are some systematic effects included which are not yet quantified. In particular, 
there are apparent edge effects from where the apertures extend beyond the observed field 
which produce spurious structures at the border of the observed field. These edge effects 
seem to also skew some of the observed mass peaks towards the field edge. However the 
persistence of structures despite varying apertures and polynomial indices, as well as their 
correlation with known mass structures (such as the visible subclustering of cluster member 
galaxies) indicate that these structures are sourced by actual physical structures. In the 
context of these substructure reconstructions, a mass-sheet transformation k — > Xn + 1 — A 
is equivalent to K — > K + In |A|, or simply a constant offset in K. 
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6. Intrinsic Flexion 

An important consideration for including flexion into combined strong and weak lensing 
mass reconstructions in future work is quantifying the amount of flexion-like shape intrinsic 
to unlensed galaxy profiles. Just as real galaxy profiles have a distribution of non-zero 
ellipticities which produce a scatter in the shear inferred from individual objects, there 
are also distributions of the third-order moments in unlensed galaxy profiles which create 
an "intrinsic flexion" signal. This unlensed-shape-induced scatter is best characterized by 
the scatter in the dimensionless products a bs^i and a b s ^3, where a b s is the observed, 
image plane size of the image. The intrinsic flexion scatter should not be correlated with 
the astrophysically-sourced lensing flexion, so each component of the flexion fields can be 
used as an independent probe of the scatter along a single coordinate axis. 

Table [2] lists the mean and standard deviation of these dimensionless products a b s $ mn 
for the flexion fields measured in A1689, characterized in three different object samples for 
each flexion field. The full sample includes all objects used for the flexion analysis in §[5| 
"Low error" and "high error" samples are defined independently for 1-flexion and 3-flexion 
based on the sum in quadrature of their respective flexion parameter error estimates: the 
full sample is divided about the median 1-flexion error of 0.029 arcsec -1 and the median 
3-flexion error of 0.049 arcsec -1 . The scatter we observe is similar to values previously 
reported, though caution is warranted when directly comparing numerical values, due to 
the difference between our treatment of the mass-sheet degeneracy and that of Goldberg fc 



Leonard (2007) 



Another notable feature is the higher 1-flexion scatter in the low error subsample than 
in the high error sample, which we attribute to the astrophysical flexion signal. The low 
error sample are those objects most likely to include a non-zero measured flexion value. The 
objects which have larger true flexion field values are more likely to be well-constrained by 
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the AIM fitting, and thus will more likely be in the low error sample and increase the scatter 
in our unitless figure-of-merit. Sub-dividing the full sample by their radius from the cluster 
center has no discernible effect on the scatter, which is expected. The flexion signal is less 
dependent on the overall cluster potential, and more dependent on the local substructure. 
A detailed study of the scatter in the measured flexion fields and the relative contributions 
from galaxy shape, galaxy-galaxy lensing, and cluster substructure lensing, is beyond the 
scope of this paper, but will be addressed in future work. 3-flexion measurements are 
typically noisier than 1-flexion, and the increase in scatter in the high noise sample is an 
indication that the 3-flexion scatter may be dominated by the measurement noise. 

7. Summary and Discussion 

Obtaining detailed measurements of the mass structure in galaxy clusters is an 
important step in understanding these important cosmological objects. To this end, a 
variety of mass reconstruction techniques based on measuring the lensing distortion of 
background galaxies have been developed. The combination of multiple orders of lensing 
effects, such as weak lensing shear with strong lensing measurements, allows for a mass 
reconstruction over a large area and a range of lensing field strengths. 

However, near the center of the cluster the convergence can vary significantly over the 
scale on which galaxy ellipticities are typically averaged to extract the weak-lensing shear, 
increasing the scatter in the shear measurements and reducing the resolution of the weak 
lensing mass maps in these areas. It is not, however, possible to measure a strong lensing 
signal for all background galaxies in this regime since not all objects will be multiply 
imaged. An alternative approach is necessary to extract the lensing information from these 
distorted background images. Second-order weak lensing distortions, called flexion, are high 
signal-to-noise probes of mass substructure on small angular scales. Flexion measurements 
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alone can be used to measure mass structure, and including them in joint weak+strong 
lensing measurements of cluster mass distributions will better resolve cluster substructures. 

In recent years, two main techniques for measuring flexion (shapelets and HOLICs) 
have been applied to measure mass structure in the galaxy cluster Abell 1689. We have 
developed an alternative method, the Analytic Image Model (AIM) method. This approach 
has several desirable properties for flexion measurement. Because the AIM method is based 
on model optimization, a full covariance matrix including both lensing and non-lensing 
parameters accompanies the flexion field estimates for each observed background galaxy. 
This provides direct error estimates for all components of the flexion fields as well as an 
indicator of parameter degeneracies. Additionally, the AIM method is fully invariant to 
the mass-sheet degeneracy, a property that previous flexion analyses of observational data 
have not had. And because the AIM method forward-models quadratic lensing rather 
than inverting measured properties to obtain flexion estimates, it provides an important 
systematic check on the shaplet and HOLIC based methods. 

Using simulated data images across a broad span of parameter space, we have 
demonstrated that the AIM method accurately reproduces both the lensing parameters 
and the intrinsic shape parameters. By fixing the shear parameters to an estimate of the 
true shear field, we recover the flexion parameters. We also recover the intrinsic shape 
parameters up to the error in the ellipticity induced by an erroneous fixed shear value. 
These simulations validate the accuracy of the AIM method for measuring flexion. 

We applied the AIM method to measure flexion fields in the galaxy cluster Abell 1689, 



and used a modified version of the 1-flexion aperture mass statistic presented by Leonard 



et al. (2010), accounting for the mass-sheet degeneracy, to infer mass structures which are 
consistent with previous measurements. The structures we observe are significantly detected 
using multiple aperture functions, meaning that they are sourced by physical structures in 
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the cluster. This application both highlights the utility of flexion measurements for mass 
reconstructions in clusters, as well as the validity of the flexion measurements obtained via 
the AIM method. 

A robust literature of weak lensing shear and strong lensing mass measurements already 
exists, and many studies combine both weak and strong lensing to improve constraints on 
the mass distribution in clusters of galaxies. Including flexion in these multi-scale mass 
reconstructions will allow for a more detailed mapping of the mass distribution in the 
central parts of clusters, regions where the linear weak lensing approximation does not 
capture all of the measurable lensing information and strong lensing data is not ubiquitous. 
A complete comparison between the AIM method and the two previously described methods 
remains to be done; such a comparison will be able to characterize the noise properties 
and systematics of each technique, as well as their relative strengths and weaknesses. 
Characterization of the noise from flexion-like shape properties in unlensed galaxy images 
and an understanding of how that noise contaminates the lensing flexion signal is needed 
for accurate mass reconstructions. Due to its fundamentally different approach, the AIM 
method will provide an important constraint on the intrinsic noise of flexion measurements, 
separating out the systematic effects introduced by measurement techniques. 
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2834-MIT-SAO-4018 issued by the Chandra X-Ray Observatory Center, which is operated 
by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract 
NAS8-03060, and through grant HST-GO-11099 from the Space Telescope Science Institute 
(STScI), which is operated by AURA, Inc., under NASA contract NAS5-26555 and 
NNX08AD79G. 
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Fig. 1. — A comparison of fit and input flexion values from simulations for 1-flexion and 
3-flexion. We plot only one component of each flexion field, but the other is very similar. 
The points in both panels are drawn from the same representative subset of the simulated 
objects, and the overlaid red line is a unity line, not a fit. As is typical, the 3-flexion errors 
are larger and the fit values are more scattered. The angular flexion units assume the HST 
ACS pixel scale. 
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Table 1. Allowed ranges for model fit parameters. L is the full side length of the 
single-galaxy data image being fit, which is chosen to be 1.5 times the observed semi-major 
axis length. Angular units assume a HST ACS WFC pixel scale of CC.'05/pixel, as in §[5] 



Parameter 


Significance 


Range 




Reduced 1-flexion 


—2 to 2 arcsec -1 


*31, *32 


Reduced 3-flexion 


—2 to 2 arcsec -1 


log S 


Total flux 


-10 to 10 


a 


Image Size 


0.05 to 3L arcsec 




Image Center 


—L to L arcsec 




Ellipticity 


-1 to 1 



Table 2. Mean and standard deviations for the products using the SExtractor-determined 
image sizes in the full, low error, and high error object samples. 



Sample (a* ln ) a a<E , ln (a^ 3n ) a a ^. in 

Full 0.000 0.042 0.000 0.048 

Low Error 0.004 0.047 0.000 0.041 

High Error -0.004 0.037 -0.001 0.055 
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Fig. 2. — Example A1689 fit objects. In all three rows, the left panel is the A1689 F775W 
data image used for the image analysis, the middle panel is the best-fit model image, and 
the right panel is the residual image. The top row is a well-fit object with y 2 /~D.o.¥ .< 1.5, 
the middle row is an object just below the threshold x 2 /D-0-F.=1.5, and the bottom row is 
an object with x 2 /D.o.F.> 1.5, and thus rejected from the mass reconstruction. 
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Fig. 3. — Aperture weight (left) and 1-flexion kernel (right) as a function of scaled radius 
for the polynomial indices used in the K ap reconstruction. Q and r have been scaled by 
the aperture radius R to make them unitless quantities. A higher polynomial index implies 
sensitivity to a narrower range of radii within the aperture. The kernel is negative, meaning 
that 1-flexion measurements oriented towards the aperture center contribute positively to 
the aperture statistic. 
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Fig. 4. — K ap signal-to-noise ratio (SNR) contours for R = 60" and three polynomial indices. 
Top left: I = 3; top right: I = 5; bottom left: I = 7; bottom right: four-filter coadded image 
without contours. The SNR contours beginn at SNR=5 and increase in steps of 2. Labels 
indicate celestial north/east and the angular scale. 1"=3.1 kpc at the cluster redshift. 
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